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Abstract 

The system of particles interacting via multibody interatomic potential of 
general form is considered. Possible variants of partition of the total force acting 
on a single particle into pair contributions are discussed. Two definitions for 
the force acting between a pair of particles are compared. The forces coincide 
only if the particles interact via pair or embedded-atom potentials. However in 
literature both definitions are used in order to determine Cauchy stress tensor. 
A simplest example of the linear pure shear of perfect square lattice is analyzed. 
It is shown that, Hardy's definition for the stress tensor gives different results 
depending on the radius of localization function. The differences strongly depend 
on the way of the force definition. 
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1 Introduction 

Classical molecular dynamics (MD) [U12], based on the numerical solution of the New- 
tonian equations of motion of many interacting particles, has been widely used for 
physical modeling and simulation for several decades. Interactions between particles 
are usually described by so-called empirical interatomic potentials. A variety of po- 
tentials has been proposed: starting with simple pair potentials |3j, more accurate 
embedded-atom potential [I] and ending with complex bond-order potentials [5J [BJ [Z] • 
The MD simulation technique requires calculation of the total force acting on every 
particle. If the potential energy is known then force Fj acting on particle number i is 
determined by the following formula [8] 

^({R,}f =1 ) 

F * = dR t ' (1) 

where {Ri}^ 1 is the set of radius- vectors of particles; A" is the total number of particles; 
U is the total potential energy of the system. Thus formula ([1]) is sufficient for MD 
simulation. However the problem of interpretation and verification of results of the 
simulation is not so straightforward. The usual way is to calculate continuum variables, 
such as stress tensor and heat flux, during MD simulation and compare the results with 
the predictions of continuum theory. This problem was addressed by many authors 
starting from the late 19th century [9]. Comprehensive reviews on this topic may 
be found in papers [TUl PIT] . The majority of approaches use different energy and 
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forces assumptions. Though the correctness of these assumptions is clear for pair 
potentials, for multibody potentials, such as Tersoff potential, this is no longer the 
case. In particular, in Hardy's formalism [12] it was assumed that the total force can 
be expressed by the summation of pair contributions satisfying Newton's third law, 
i.e. 

Fj = Fjj , Fjj = — Fji. (2) 

Here and below the summation is carried out over all particles in the system. This 
assumption was examined in paper [13]. It was shown that it holds for systems with 
three-body forces [5JEJE]. The generalization for N-body forces was carried out in [T3] , 
In paper [llj it was stated that the partition ([2]) can always be carried out. However, 
the physical meaning of F^ is clear only in the case of pair potentials. One can easily 
understand the reason for this ambiguity. Systems with pair interactions are simply 
equivalent to sets of particles connected by longitudinal springs. Therefore F^ is just 
a force caused by the deformation of the spring. In contrast the simplest systems with 
multibody interactions can be imagined as a set of particles connected by longitudinal 
and angular springs between the bonds. In this case every angular spring belongs to 
three particles and causes the forces acting on all of them. Even in this simple case 
partition (j2J) is not straightforward. Note that commonly used multibody potentials, 
like Tersoff [6], are even more complex. Another definition for Fy was given in pa- 
pers [15] and [16]. It was stated in \Tjg\ that division (J2J may be non-unique. However 
this statement was not proved. 

In the present paper the problem of the partition of the total force into pair contri- 
butions is discussed in detail. Two definitions for the force Fy acting between a pair of 
particles (material pointsj]] are considered. It is shown that the partition is not unique 
and does depend on the way the potential energy is represented. The influence of the 
partition on the value of stresses is analyzed for the simplest example — a linear pure 
shear of a square lattice. 



2 Comparison of different definitions for the force 
acting between two particles 

Let us begin with consideration of the results obtained in papers [T3|fH]. The problem 
of calculation of in the case of three-body forces was addressed in [13]. It was shown 
that both the Stillinger- Weber [5] and the Tersoff [6] potentials can be expressed in the 
following form 

U — 77 Uj, Ui = Uijk{Rij, Rik, Rkj)i Rij — (3) 

i j^i-kjti,j 

Here and below = Rj — Rj. Then according to definition ([T]) the total force acting 
on the i-th particle is 

Tp _ \ ^ -rp rp de f _ 1 \ ^ djJJjjk + Ukij + Ukji + Ujik + Ujkj + Ujki) / , x 



1 Note that particles with rotational degrees of freedom, internal structure, etc. are not considered. 
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Here and below == Rij/Rij. The definition for force Fjj acting between particles 
i and j arises in a natural way while calculating the derivative. On the other hand 
in paper [13J it was shown that the same results can be obtained using the following 
definition 

def BU 

" - Oil,/'"' [ } 

One can see that equations ([2]) are satisfied and furthermore are central, i.e. par- 
allel to the vector connecting particles i and j. The above mentioned approach was 
generalized in paper [14] . Let us assume that the total potential energy depends on all 
the interatomic distances in the system, i.e. 

U = U{{R kn } k>n>k ). (6) 

Formula ([6]) is the most general form for potential energy of the atomic system. Ac- 
cording to definition ([1]) the total force Fj is 

F. = -T—^«=-T—e~ (7) 
2^ dR dR . ^dR^- [ ) 

k,j>k J j^i J 

One can see that again the definition (JS]) naturally follows from the derivation. Thus 
in the general case ([6]) the total force can be expressed in form (J2J). 

The partition discussed above requires the potential energy U to be represented 
as a function of all the interatomic distances in the system (jH]). The requirement 
is automatically satisfied in the case of pair potentials [5] and the embedded-atom 
potential [4]. However representation fl6]) may be inconvenient in the case of bond- 
order potentials, such as Stillinger- Weber j5] and Tersoff [6], which depend on the 
angles between bonds. Let us consider the different approach proposed in papers [T5] 
and [16J. Assume that the total potential energy of the system has the form 



Obviously all potentials mentioned above can be represented in the form ([H]). In con- 
trast to the previous approach, the geometry of the atomic system should be represented 
via vectors {Rjjjj^j, but not interatomic distances. At first glance it seems that both 
approaches are equivalent. Let us show that in general it is not true. Calculating the 
total force Fj using the definition (TT]) and expression (JSJ) one obtains 

. \ - \ - 9Uj_ dRjk = y> ( dUj_ _ dUA d= ef 9Uj_ _ dU^ 

j tit 1 ] jt 11 

One can see that introduced definition for Fy satisfies Newton's third law, i.e. Fjj = 
—Fjj. Therefore F^- may be considered as a force acting between particles i and j. 
It can be shown that, in the case of pair potentials and embedded-atom potentials, 
forces F^ are central. Moreover in these particular cases definitions ([5]) and exactly 
coincide. However in general this is not true. The first distinction and, probably, the 
most important one is that the partition of the total energy ([S]), used in the expres- 
sion (jnD, is not unique. According to formula §§§ the partition can, in general, affect the 
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value of the force Fy. However the majority of commonly used potentials are based on 
partition (jHJ) as well. For example, in the case of Tersoff and Stillinger- Weber potentials 
the partition is determined by formula ([3]). The comparison of different partitions will 
be carried out below for a simple example. 

Let us compare definitions ([SD, (JSJ) in the case of three-body potentials. Rewriting 
the expression for potential energy ([3]) in the form analogous to (jHJ) one obtains 



U ^ Uj, Ui ^ Uikn (Rjfc, Rjn) , ?7j/ cn (Rjfc, Rjri) Uikn(Riki Rim Rkn) 



k,n^k 



One can show that in this case the expression ([9]) for F»j takes the form 



dcf 



2 ^ 



d 



dUji 







dUij 



1 1- j 



(10) 



(11) 



Let us substitute the last formula from (flOj) into (TlT]) and calculate the derivatives. 
Then Fij can be expressed in the form analogous to (j4j) 



d (Uijk + Uikj + Ujik + U. 



jki, 



dRi-; 



d (Ujik + Ujki) d (U^k + Uikj) 

'^ik i" oi T-> ^"kj 



dRik 



OR 



(12) 

One can clearly see that expressions @ and ffT2l are different. In particular, from (jlj) 
it follows that forces are parallel to Rjj. According to ffl2|) this is not true. Further 
these statements will be explicitly shown for a simple example (see equations f[T4"|) . 

dm (Hi). 



3 Calculation of Cauchy stress tensor. The sim- 
plest example 

Let us consider the practical consequences of different ways of partitioning of the to- 
tal force. It was mentioned above that results of MD simulation are usually compared 
with predictions of continuum theory. Equivalent continuum quantities, such as Cauchy 
stress tensor, are calculated for this purpose. The stress tensor characterizes external 
forces acting on the material surface surrounding the volume selected from the contin- 
uum media. Therefore the forces should be divided into internal and external in order 
to calculate stress. The analogous situation occurs if one defines equivalent stress tensor 
for discrete system. Obviously it is impossible to express stress tensor in terms of the 
total forces acting on the particles. Therefore partition ([2]) is required. It was shown 
above that the partition is not unique. Thus we obtain the result that the stresses in 
a material can depend on the choice of definition for F^. Let us consider the simplest 
example, notably linear pure shear of a perfect zero-temperature square lattice. Let 
particles interact via angular springs connecting two neighboring bonds. Only nearest 
neighbors will be taken into account. Longitudinal springs may also be added but their 
contribution to forces F^ is the same for both expressions 0. Let us consider 
particle number with radius- vector R and denote its neighbors as in figure 1. Note 
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Figure 1: Particle number 
and its nearest neighbors. 



that this numbering takes symmetry into account, for example, Rqi = — Ro(-i)- Let 
us assume that potential energy of the spring connecting bonds Roi and R02 is given 
by the following formula 

U m = °- (e 01 • e 02 ) 2 = °- ( R oi+fo2-R 2 n \ \ (13) 



2 2 \ 2R Q1 R 02 / 

In the case of small deformations formula ffT3]) corresponds to the energy of harmonic 
angular spring with stiffness c. Let us derive expressions for forces F i, F 2, Fi 2 , 
Fx(_ 2 ). The remaining forces can be obtained using symmetry or Newton's third law, 
for example, F 12 = F(_ 2 )(-i)> F i = — F 10 . First let us use definition (jSJ). Substituting 
formula (fT3"l) into formula (jl]) and linearizing the resulting expression in the case of 

small shear (tp = arcsin(e i • e 02 ), y<l) one obtains 

^ 2^2cp 2y/2ap . , 

Fqi^O, F 02 ~0, F 12 ^^^e 12 , F 1( _ 2) « -^^e 1{ _ 2) . (14) 

where R = i?oi = Ro2- Here and below terms of order higher then ip are neglected. It 
is clear that forces (1141) are central. 

Let us use definition (JHJ). It was mentioned above that representation (JSJ) of the 
total energy used in formula ([9]) is not unique. In the case under consideration particle 
number is surrounded by four nearest neighbors and four corresponding spring. Let us 
compare two possible ways of the partition of the total energy. First, let us assume that 
the energy of the springs contribute to energy Uq of particle number only. This way of 
partition is the simplest one and it is similar to the partition used in the definition for 
Tersoff and Stillinger- Weber potentials ( TTUj) . Calculating forces F i, F 02 , Fi 2 , Fi(_2) 
in the framework of the assumptions used during the derivation of formula (lT4"j) one 
obtains 

-^-e 02 , F 02 ~ — — 

Obviously the forces (j!5p are non-central and differ from the previous result ( I14p . Note 
that in formula f JT5|) forces Fi 2 ,Fi(_ 2 ) are exactly equal to zero. The second way of 
partition of the total energy was proposed in paper [16]. Energy t/012 of the spring was 
divided between particles number 0, 1 and 2 in equal portions. Calculating the forces 



01 ~ B-e 02 , F 02 w — e i, F12 = Fi(_ 2 ) = 0. (15) 
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one obtains 



8 op 



3R 



©02, Fq2 



8ctp 



3R 



eoi, F12 



3R 



ei2, Fi(_ 2 ) 



3R 



ei(_ 2 ). 



(16) 



Note that in contrast to formula (IT5|) . forces F 12 ,Fi(_ 2 ) determined by formula f fTBj) 
are not equal to zero. 

Let us calculate shear stress acting on cross-section with normal eoi of the crystal. 
According to classical continuum mechanics definition, the shear stress acting on the 
cross-section is equal to the component of the force parallel to vector e 2 per unit 
length. One can prove using formulas (THj) . (TT51) . (TT6"j) that in all the cases the absolute 
value of shear stress r can be expressed as 



The second formula from (1TTI) coincides with the well-known expression for elastic 
constant C44 of the square lattice, obtained, for example, in [17] . Index star is used 
in order to mark the exact solution. Thus, at least in the case under consideration, 
stresses calculated using different expressions for F^- coincide. 

Application of mentioned above definition for the stress is not very convenient, 
especially for dynamical problems involving large thermal motion, structural transfor- 
mations, fracture, etc. Therefore in practice different approaches for calculation of 
stress are used (see papers [TUl [EE] for detailed reviews). In particular, Hardy's for- 
malism [T2] is frequently used [TTj [T5| ITE] . In papers [T31 [EH] it was shown that in the 
framework of Hardy's formalism the potential part of Cauchy stress tensor at spatial 
point Rq has the form 



where ip is so-called localization function (see [15] for details). Formula (TTH|) was used 
in paper [13J and in the Appendix of paper [16J. However the quantity Fij has different 
meanings in them. Formula ([5]) was used as a definition for Fy in paper [13]. In 
contrast in paper [16] quantity F^ was calculated using formula ([9]). It was shown 
above that, in general, forces calculated with the use of (EJ) and are different. Thus 
corresponding stress tensors are also different. The accurate comparison of the stress 
tensors will be addressed in a separate paper. In the present work only the simplest 
example is analyzed. 

Let us return to the example of linear pure shear of the square lattice. Consider 
Cauchy stress tensor at the point Ro, where particle number is placed. One can prove 
substituting formulas ( TH| . ( TH>j) . ( TTB"|) into formula ( TTHj) that in all these cases stress 
tensor is symmetric. Note that, in general, this is not obvious as forces (j!5p and (1161) 
are non-central. Let us calculate elastic constant C44 of the system as it completely 
determines stresses in the case under consideration. The exact solution of this problem 
is given by the second formula from ffTTj) . For simplicity the radial step function was 
used as a localization function, i.e. ip = l/(7rr|) inside the localization volume and 
zero in the remaining space, where is radius of the localization volume. The elastic 



t = C^\(p\, C\ 



.* def 4c 
44 - ^2- 



(17) 




(18) 
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Figure 2: Calculated elastic constant C 44 divided by the exact solu- 
tion Cl A . Squares, triangles and circles correspond to elastic constant 
calculated using expressions (CHI) . f|T5|) . (fTBj) respectively. 

constants calculated with different values of are shown in figure 2. One can see from 
figure 2 that for small r^/R elastic constants calculated using expressions (J5J) and (jH]) 
for Fjj are different. Furthermore in the second case elastic constant C 44 depend on the 
partition of the total energy. However the expressions converge to the same value C| 4 
with increasing radius of localization function. Note that the elastic constant which 
corresponds to definition ([5} converges to the exact solution more rapidly then the 
elastic constant which corresponds to definition (Q. The practical consequence of this 
fact is that the first approach requires a smaller value of Vl than the second one and 
it is therefore more efficient from the computational point of view. 

4 Results and discussion 

Let us summarize the results. The problem of the partition of forces acting in the 
discrete system with multibody interactions into pair contributions is analyzed. Two 
methods for partitioning, differing in the way representation of potential energy U, are 
discussed. It is shown that in the framework of both methods the total force Fj is 
expressed as a sum of pair contributions Fy. In both cases the definition for F^ is 
rather natural as the analog of Newton's third law is satisfied, i.e. F^ = — Fjj. However 
the definitions coincide only in the case of pair potentials and the embedded-atom 
potential. In particular, according to the approach proposed in papers [131 EH]) f° rces 
Fij are always central. In the framework of the approach proposed in papers [T5| [16] 
forces are, in general, non-central. The difference was explicitly shown in the case of 
three-body potentials. Thus the partition mentioned above is not unique. It depends on 
the representation of the potential energy. The influence of the partition on the value 
of stress tensor is analyzed for the simplest example of linear pure shear of perfect 
zero-temperature square lattice. Two definitions for stress are considered. If stresses 



are calculated as a force acting on the cross-section of the crystal per unit length, all 
expressions for F^- considered above lead to exactly the same value of stresses. In 
contrast, Hardy's definition (JTHjl gives different values of stresses depending on the 
size of localization volume and the differences depend strongly on the way of force 
definition. Accurate comparison of different expressions for stress tensor in the case on 
nonlinear deformation involving thermal motion is addressed in our future work. 

Finally let us note that in the case of multibody interactions a common point 
of view on the definition for the force acting between two particles does not exist. 
Summarizing the results of the present paper one can formulate several advantages of 
definition (jSJ). First, definition (jSj) does not explicitly depend on partition (JSJ) of the 
total energy, which is, in general, not unique. Secondly, it was mentioned that forces 
defined by formula d^J) are central. In this case equation of moment of momentum bal- 
ance is satisfied for any subsystem of the discrete system [18J. It leads, in particular, 
to unconditional symmetry of corresponding Cauchy stress tensor. Note that in the 
framework of classical continuum mechanics the symmetry is required [18J. On the 
other hand, in paper [16] it was shown that definition is more appropriate for the 
formulation of equivalent micromorphic continuum theory for discrete systems. The 
third appealing feature of definition © is that, at least in the considered case, the elas- 
tic constant C 44 and therefore the stress tensor, calculated with the use of formula (jSJ), 
converges to the exact solution more rapidly than in the case of definition (jUJ). In the 
last case it is shown that the convergence rate depend on the way of partition of the 
total energy between particles. 

The author is deeply grateful to Prof. Y. Chen, Prof. Wm.G. Hoover, Prof. E.A. 
Ivanova, Prof. A.M. Krivtsov and Prof. J. A. Zimmerman for useful discussions and 
inspiration. 
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